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ABSTRACT 


The  periodogram,  the  square  of  the  magnitude  of  the  Fourier 
Transform,  is  widely  used  to  estimate  the  spectral  content  of 
sampled  processes.  The  performance  of  the  periodogram  is 
degraded  by  spectral  leakage.  This  is  the  consequence  of 
processing  finite-length  data  records.  Classical  means  of 
enhancing  periodogram  performance  are  the  use  of  tapered 
window  functions  and  averaging  of  several  periodograms .  These 
methods  smooth  the  spectral  estimate,  but  at  a  loss  of 
resolution.  A  non-stationary  Kalman  filter  was  applied  to  the 
periodogram  of  untapered  (i.e.,  rectangular  windowed)  time 
data  in  an  effort  to  smooth  the  noise  portions  of  the 
periodogram  while  leaving  the  main  spectral  response 
unaltered.  The  Kalman  filter  was  able  to  enhance  the 
periodogram.  Best  results  were  obtained  in  the  single 
spectral  peak  case.  Even  in  the  case  of  multiple  spectral 
peaks,  the  resolution  of  the  unfiltered  periodogram  was 
largely  preserved  since  the  filtering  algorithm  was  designed 
to  selectively  smooth  the  noise-only  segments  of  the  spectral 
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I .  INTRODUCTION 


The  periodogram,  the  square  of  the  magnitude  of  the 
Fourier  Transform,  is  widely  used  to  estimate  the  spectral 
content  of  sampled  processes.  The  periodogram  remains  popular 
in  the  face  of  more  modern  spectral  estimation  techniques 
(i.e.,  parametric  modeling)  due  to  its  low  cost  and  ease  of 
implementation  in  real  time.  The  performance  (ability  to 
detect  signals  in  noise)  of  the  periodogram  is  degraded  by 
window  function  sidelobe  effects.  This  is  the  unavoidable 
consequence  of  processing  data  records  of  finite  length.  In 
addition,  the  periodogram  may  have  a  fairly  large  variance 
(i.e.,  mean  equals  the  standard  deviation  under  noise-only 
conditions) .  A  classical  means  of  enhancing  the  performance 
of  the  periodogram  is  the  use  of  tapered  window  functions, 
such  as  the  Hamming  window,  in  order  to  minimize  the  effects 
of  the  discontinuity  at  the  boundaries  of  the  finite 
observation.  Another  common  method  is  to  average  a  series  of 
periodograms  in  an  effort  to  smooth  the  spectral  estimate 
(i.e.,  reduce  the  variance  of  the  estimate).  Almost 
invariably,  the  consequences  of  these  techniques  are  a 
broadening  of  the  main  spectral  peaks  and  a  corresponding  loss 
of  spectral  resolution.  What  is  proposed  here  is  an 
application  of  a  non-stationary  Kalman  filter  to  the  sequence 
presented  by  the  periodogram  of  untapered  (i.e.,  rectangular 
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windowed)  time  data.  The  objective  is  to  filter  (smooth)  the 
noise  portions  of  the  spectral  estimate  and  leave  the  main 
spectral  responses  unaltered.  The  result  is  that  the  dominant 
spectral  peaks  will  be  highlighted  against  the  noise  "floor" 
out  of  which  they  rise.  Since  the  main  spectral  peaks  are 
unaltered,  the  resolution  of  the  original  periodogram  is 
preserved.  Using  the  test  cases  of  single  and  multiple 
sinusoids  in  Gaussian  white  noise,  the  Kalman  filter's 
performance  was  evaluated  for  signal  detectability  and 
resolution  at  different  input  signal-to-noise  ratios  on 
multiple  noise  realizations.  The  effects  of  varying  the 
filter's  detection  parameter  and  the  data/transform  length 
were  also  investigated. 
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II.  CLASSICAL  SPECTRAL  ESTIMATION 


A.  BACKGROUND 

Estimation  of  the  power  spectral  density  (PSD)  of  sampled 
deterministic  or  stochastic  processes  is  usually  based  on 
techniques  employing  the  Fast  Fourier  Transform  (FFT) .  These 
techniques  are  computationally  efficient  and  produce  good 
results  for  many  different  types  of  signals.  There  are, 
however,  two  significant  limitations  associated  with  the  FFT- 
based  techniques.  First  and  foremost  is  the  problem  of 
frequency  resolution,  that  is,  the  ability  to  distinguish 
between  the  presence  of  one  or  several  spectral  components  in 
a  given  sample  set  of  data.  Frequency  resolution  of 
stationary  signals  varies  with  the  specific  technique  employed 
but,  in  general,  it  is  proportional  to  the  reciprocal  of  the 
time  interval  represented  by  the  sample.  The  second 
limitation  of  the  FFT-based  methods  is  caused  by  the  windowing 
of  the  data  that  occurs  during  processing.  Windowing  causes 
"leakage"  in  the  spectral  domain.  Energy  in  the  main  lobe  of 
a  spectral  response  "leaks"  into  adjacent  sidelobes,  obscuring 
and  distorting  the  spectral  responses  due  to  other  frequency 
components  that  may  be  present.  In  some  cases,  weak  spectral 
responses  may  be  completely  masked  by  the  sidelobes  of 
stronger  spectral  responses  and  thus  go  undetected.  Careful 
selection  and  use  of  tapered  data  windows  can  reduce  sidelobe 
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leakage,  but  always  at  the  cost  of  reduced  frequency 
resolution  [Ref.  1] . 

B.  CLASSICAL  SPECTRAL  ESTIMATION  TECHNIQUES 

The  two  best-known  classical  spectral  estimation 
techniques  are  the  Blackman-Tukey  method  and  the  periodogram. 
The  Blackman-Tukey  approach,  introduced  in  1958  [Ref.  2], 
first  estimates  the  autocorrelation  function  from  the  data 
and  then  Fourier  transforms  the  correlation  estimates  to 
obtain  a  power  spectral  density  estimate.  The  Blackman-Tukey 
spectral  estimator  is  given  by: 

N-\ 

/07  (/)  =  2n/k)  (2.i) 

where 

,  N-I -* 

—  rx-(n)x(,nk);  k  =  0, 1,2,...,(N-1) 

U*)“  N  »»o 

/«(  -*);  1HN-IHN-2), . -1  .  (2.2) 

This  is  a  biased  estimator  of  the  true  autocorrelation 
function  since: 

r(U  ,  ,2.3) 

The  mean  value  of  the  autocorrelation  function  estimator 
shows  that  a  triangular  (Bartlett)  window  is  applied  to  the 
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true  autocorrelation  function.  It  is  possible  to  use  an 
unbiased  autocorrelation  function  estimator  by  replacing  the 
normalization  by  1/N  in  (2.2)  with  l/(N-|k|).  This,  however, 
can  lead  to  a  negative  spectral  estimate  since  the  unbiased 
autocorrelation  estimator  does  not  guarantee  a  positive  semi- 
definite  sequence.  The  Blackman-Tukey  approach  was  the  most 
popular  spectral  estimation  technique  until  the  introduction 
of  the  FFT  algorithm  [Refs.  3  and  4). 

The  periodogram  spectral  estimate  is  obtained  from  the 
square  of  the  magnitude  of  the  Fourier  transform  of  the  data. 
The  data  may  be  weighted  by  a  window  function  and/or  zero- 
padded.  The  true  spectral  estimator  is  given  by: 


'«(/)=  x 


1 


2M  +  1 


A1 


I>(»)exp(-/2/r/>i)| 


(2.4) 


If  we  ignore  the  expectation  operator  and  use  only  the 
available  data,  the  spectral  estimator,  denoted  as  the 
periodogram,  is  given  by: 


^  V'FK  (/  )  - 


J 

~N 


£.t(?j)exp(-/2/r/») 

rt- 0  • 


(2.5) 


The  periodogram  produces  best  results  when  an  integer 
multiple  of  periods  of  constituent  frequency  components  arc 
present  in  the  observation.  Despite  the  advent  of  more  modern 
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techniques,  the  periodogram  remains  a  popular  means  of 
spectral  estimation  because  it  can  be  easily  and  inexpensively 
implemented  in  real  time. 

In  general,  the  Blackman-Tukey  and  the  periodogram 
spectral  estimates  are  not  identical.  If,  however,  the  biased 
autocorrelation  estimate  (2.2)  is  used  and  as  many 
autocorrelation  lags  as  data  samples  (N)  are  computed,  then 
the  Blackman-Tukey  and  periodogram  estimators  yield  identical 
numerical  results. 

C.  WINDOW  FUNCTIONS 

Every  set  of  data  is  finite  in  duration.  Processing  a 
finite  duration  observation  presents  special  problems  to  the 
harmonic  analysis  of  the  data.  Some  considerations  should  be 
given  to  detectability  of  spectral  components  in  the  presence 
of  nearby  strong  components  and  their  resolvability.  Let  the 
data  to  be  processed  consist  of  N  uniformly-spaced  samples  of 
the  observed  signal.  The  FFT,  the  basis  of  the  periodogram 
spectral  estimator,  assumes  sequences  to  be  periodic.  In 
other  words,  the  sample  set  under  analysis  is  assumed  to  be 
one  complete  period  of  an  infinitely  long  periodic  sequence. 
The  selection  of  a  finite  time  interval  of  NT  seconds,  where 
T  is  the  time  between  samples,  and  of  the  orthogonal 
trigonometric  basis  over  this  interval  leads  to  an  interesting 
peculiarity  of  the  spectral  expansion.  From  the  continuum  of 
possible  frequencies,  only  those  which  coincide  with  the  basis 
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functions  (the  bin  centers  of  the  FFT)  will  project  onto  a 
single  basis  vector.  All  other  frequencies  will  exhibit  non¬ 
zero  projections  on  the  entire  basis  set.  This  phenomena  is 
called  spectral  leakage  and  is  a  consequence  of  processing 
finite  duration  data  records  [Ref.  1], 

Spectral  components  with  frequencies  other  than  those 
corresponding  to  the  FFT  bin  centers  will  typically  be  present 
in  the  observed  data.  Components  with  frequencies  not  at  bin 
centers  are  not  periodic  in  the  observation  window.  The 
periodic  extension  of  a  signal  which  does  not  coincide  with 
the  natural  periods  of  its  constituent  frequency  components 
exhibits  discontinuities  at  the  boundaries  of  the  observation. 
These  discontinuities  are  responsible  for  spectral 
contributions  (leakage)  over  the  entire  range  of  the  FFT 
frequency  bins. 

Since  we  are  constrained  to  deal  only  with  finite-length 
data,  we  are  forced  to  make  certain  assumptions  about  the  data 
outside  of  the  observation  interval.  The  finite  data  record 
may  be  considered  as  having  been  obtained  by  multiplying  an 
infinite  length  data  sequence  with  a  simple  rectangular 
function: 


)l  =  0,1,2,...,(N  -  1) 


otherwise 


(2.6) 
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The  assumption  that  the  data  outside  of  the  observation 
window  is  zero  is  unrealistic  but  unavoidable.  Thus,  data 
taken  "as  is"  is  actually  rectangularly  windowed.  Non- 
rectangular  window  functions  are  weighting  functions  applied 
to  the  received  data  in  order  to  reduce  the  spectral  leakage 
associated  with  finite  observation  intervals.  The  purpose  of 
the  window  is  to  reduce  the  magnitude  of  the  discontinuity  at 
the  boundaries  of  the  periodic  extension.  The  goal  of 
windowing  is,  therefore,  to  smoothly  taper  the  data  record  at 
the  boundaries. 

By  the  Convolution  Theorem,  multiplication  of  the  time 
series  by  a  window  function  corresponds  in  the  frequency 
domain  to  the  convolution  of  the  transforms  of  the  signal 
sequence  and  the  window  function.  If  we  are  using  a 
rectangular  window  and  attempting  to  detect  a  narrow-band 
signal,  such  as  a  sinusoid  in  noise,  and  the  sinusoidal 
frequency  is  not  at  a  bin  center,  the  convolution  will  spread 
or  smear  some  signal  power  into  adjacent  frequencies. 
Conversely,  if  the  sinusoid  is  at  a  bin  center,  then  we  will 
see  only  the  zero  crossings  of  the  window  transform,  and 
experience  no  leakage.  If  we  are  using  a  non-rectangular 
window  (i.e.,  a  Hamming  window),  the  convolution  operation 
will  smear  the  signal  power  into  adjacent  frequencies 
regardless  of  the  sinusoidal  frequency  being  at  a  bin  center 
or  not. 
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Leakage  has  an  obvious  negative  effect  on  the  detection 
and  estimation  of  sinusoidal  components.  Sidelobes  from 
adjacent  frequency  components  may  add  in  an  unpredictable 
fashion  to  the  spectral  peak  of  a  weak  signal,  thus  distorting 
the  power  estimate  of  that  signal.  In  extreme  cases,  the 
sidelobes  of  strong  frequency  components  may  completely  mask 
the  main  lobe  of  nearby  weaker  signals  [Ref.  3]. 

In  general,  the  convolution  of  the  window  transform  with 
the  signal  transform  means  that  the  main  lobe  width  of  the 
window  transform  is  the  limiting  factor  (in  terms  of  spectral 
response)  that  allows  separation  of  two  closely-spaced 
spectral  lines.  For  a  rectangular  window,  the  main  lobe  width 
between  the  3-dB  levels  of  the  resulting  digital  sine  function 
(the  FFT  of  a  rectangle  function)  is  approximately  the 
reciprocal  of  the  observation  interval  NT.  Leakage  effects 
can  be  reduced  by  the  use  of  windows  with  non-uniform 
weighting,  such  as  the  Blackman  and  Hamming  windows. 

Consider,  for  example,  the  problem  of  detecting  a 
sinusoidal  signal  embedded  in  Gaussian  white  noise.  Assuming 
that  the  observation  interval  does  not  contain  an  integer 
multiple  of  periods  of  the  sinusoid,  then  the  frequency  of  the 
sinusoid  is  not  at  a  bin  center  of  the  FFT.  Some  spectral 
leakage  will  occur.  Recall  from  basic  Fourier  theory  that  the 
transform  of  a  sinusoid  (say  a  cosine  function)  is  a  pair  of 
delta  functions  given  by: 

cos(2  n(S(2nf  -  2  */<,)  +  6(2nf  +  2  nf0 ))  {2'7) 
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Assuming  that  the  data  is  obtained  by  rectangular 
windowing  of  an  infinitely  long  sequence  (i.e.,  multiplication 


of  the  time  series  by  the  window  function) ,  then  the 
periodogram  will  be,  by  the  Convolution  Theorem,  the  square 
of  the  magnitude  of  the  convolution  of  the  delta  function  pair 
with  the  Discrete  Fourier  Transform  of  the  rectangle  function 
(a  digital  sine  function) .  The  digital  sine  function  is  of 
the  form: 


dn  (/)  =  T exp(-j2nfT[N  - 

sm (*fT)  . 


(2.8) 


Recall  from  Fourier  theory  that  the  convolution  of  some 
function,  call  it  F(f),  with  a  delta  function,  results  in  the 
translation  of  F(f)  to  the  location  of  the  delta  function. 
In  this  case,  the  sine  function  will  be  shifted  to  the 
location  of  the  delta  function  dictated  by  the  signal 
frequency.  If  the  location  of  the  delta  function  does  not 
exactly  coincide  with  a  bin  center  of  the  FFT,  leakage  will 
occur.  1 


At  this  point,  some  discussion  of  zero-padding  is  in 
order.  Zero-padding  the  data  sequence  prior  to  the  Fourier 
transformation  will  not  improve  the  resolution  of  the 
periodogram.  The  purpose  of  zero-padding  is  twofold.  First, 
it  will  interpolate  additional  power  spectral  density  values 
in  the  interval  [-fs/2,  fs/2],  where  f,  is  the  sampling 
frequency  [Ref.  3]  between  those  that  would  have  been  obtained 
in  a  non-zero-padded  transform.  Second,  since  the  number  of 
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observed  data  points  is  not  always  a  power  of  two,  zero¬ 
padding  is  necessary  to  make  the  sequence  length  a  power  of 
two  to  allow  the  use  of  a  FFT.  Consider  the  Discrete  Fourier 
Transform  of  an  eight-point  rectangular  window.  We  know  that 
this  transform  will  produce  a  digital  sine  function.  However, 
when  we  actually  compute  and  plot  the  transform,  we  observe 
only  a  central  spike  at  the  zero  spectral  location.  (Figure 
1) .  Why  do  we  not  see  any  of  the  side  lobe  structure  that  we 
know  must  be  present?  The  side  lobes  are  in  fact  there.  They 
are  not  visible  because  the  FFT  of  the  non-zero-padded  time 
series  interrogates  the  resultant  digital  sine  function  at  its 
zero-crossings  and  hence,  the  side  lobe  structure  is  invisible 
to  us.  In  other  words,  the  FFT  bin  centers  are  coincident 
with  the  digital  sine's  zero-crossings.  Now  examine  what 
happens  when  the  eight-point  rectangle  is  zero-padded  to 
sixteen  points  and  then  transformed  (Figure  2)  .  The  side 
lobes  are  now  clearly  visible  because  we  are  interpolating  a 
point  in  between  the  bin  centers  of  the  previous  eight-point 
(non-zero-padded)  transform.  This  principle  can  now  be 
extended  to  an  actual  spectral  estimation  example. 

Consider  a  unit  amplitude  sinusoid  embedded  in  Gaussian 
white  noise.  In  this  example,  the  number  of  data  points  N  is 
64  and  a  rectangular  window  is  used.  The  sinusoidal  frequency 
is  10.0  Hz  and  the  sampling  frequency,  f«,  is  64.0  Hz.  The 
variance  of  the  additive  Gaussian  white  noise  is  1/2000.  This 
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corresponds  to  a  signal-to-noise  ratio  (SNR)  of  30  dB  where 
SNR  is  defined  as: 


SNR  =  101og10 


(  sinusoidal  amplitude "'l 

'MY 

l  2  J 

= 101og10 

UJ 

/  noise  ) 

^variance  / 

(2.9) 


1 


where  A  =  amplitude  of  the  sinusoid. 
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MAG  OF  FFT  OF  8-PT  RECT  WINDOW 


requency  point 


In  this  example,  the  bin  centers  of  the  FFT  occur  at 
integer  multiples  of  fs/N,  which  in  this  case  is  64/64  or  1 
Hz.  Figure  3  shows  that  the  spectral  peak  is  well  defined 
since  the  sinusoidal  frequency  lies  exactly  at  a  bin  center 
and  no  zero-padding  was  performed  prior  to  transformation. 
We  do  not  see  the  side  lobe  structure  of  the  digital  sine 
(transform  of  the  rectangle  function) .  Observe  in  Figure  4 
what  occurs  when  the  frequency  detected  does  not  coincide  with 
a  bin  enter.  In  this  case,  the  frequency  is  10.7  Hz,  which 
is  clearly  not  a  bin  center.  The  side  lobes  of  the  digital 
sine  function  are  now  visible  since  we  are  not  interrogating 
the  sine  at  points  of  its  zero  crossing.  In  addition, 
spectral  leakage  has  smeared  the  signal  power  into  the 
adjacent  frequency  bins.  The  end  result  is  a  much  broader  and 
less-pronounced  main  lobe  (25  vs.  40  dB) . 

To  illustrate  the  effects  of  zero-padding,  let  us  now 
consider  the  situation  in  which  the  original  64-point  data 
record  has  been  zero-padded  to  128.  Now,  regardless  of 
whether  or  not  the  sinusoidal  frequency  is  at  a  bin  center, 
the  side  lobes  of  the  digital  sine  will  now  be  visible  as  a 
result  of  the  zero-padding  (see  Figures  5  and  6)  .  The  net 
effect  will  be  a  less  pronounced  main  lobe  due  to  the  side 
lobes.  In  the  case  of  f  =  10.7  (Figure  6),  the  main  lobe  is 
flattened  due  to  a  combination  of  the  sine  side  lobes  and 
spectral  leakage. 
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NUSOID  IN  NOISE,  BIN  CENTER 
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Figure  5.  Periodogram  of  Sinusoid  in  Gaussian  White  Noise 
f  =  10.0  Hz  (bin  center);  64  points  Zero-padded  to  128 

l  8 
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Figure  6.  Periodogram  of  Sinusoid  in  Gaussian  White  Noise; 
f  =  10.7  Hz  (not  a  bin  center);  64  points  Zero-padded  to  128 


D.  WINDOWS  WITH  NON-UNIFORM  WEIGHTING 


For  comparison,  the  original  64-point  data  records  for 
sinusoidal  frequencies  10.0  and  10.7  Hz  are  weighted  with  a 
Hamming  window  prior  to  zero-padding  and  Fourier 
transformation  (Figures  7  and  8)  .  The  Hamming  window 
function,  popular  due  to  its  good  performance  and  ease  of 
implementation,  has  a  maximum  side  lobe  level  of  -43  dB  versus 
-13  dB  for  a  rectangular  window.  The  price  paid  for  this  side 
lobe  suppression  is  increased  main  lobe  width.  The  3-dB  main 
lobe  width  becomes  1.30  bins  versus  0.89  bins  for  the 
rectangular  window.  The  Hamming  window  is  only  one  of  many 
such  functions.  An  exhaustive  comparison  of  window  functions 
and  their  use  in  spectral  analysis  is  given  by  Harris  [Ref. 
1].  Many  other  windows,  with  even  more  dramatic  reduction  of 
side  love  levels,  are  possible.  In  all  cases,  however,  the 
side  effect  is  always  a  broadening  of  the  main  lobe  with  its 
associated  reduction  in  spectral  resolution  [Refs.  1  and  3]. 


20 


gp  spnjiunnuj 


Figure  8.  Periodogram  of  Sinusoid  in  White  Gaussian  Noise; 
f  =  10.7  Hz  (not  at  bin  center);  64  points  Zero-padded  to  128; 
Hamming  window 


E.  STATISTICAL  PROPERTIES  OF  THE  PERIODOGRAM 


Consider  a  data  record  of  samples  of  Gaussian  white  noise 
having  zero  mean  and  variance  a„J.  The  periodogram  of  this 
data  will  have  a  distribution  which  is  chi-squared  with  two 
degrees  of  freedom.  The  reason  for  this  is  that  the  sampled 
Gaussian  random  process,  denoted  as  x(n) ,  has  the 
distribution: 


x{n)~N{Qfo J) 


(2 . 10) 


For  simplicity  let  us  assume  that  the  Fourier  Transform 
of  x(n)  is  normalized  by  1/SQRT(N)  ,  where  N  is  the  size  of 
the  transform.  Since  the  real  and  imaginary  parts  of  the 
Fourier  Transform  of  x(n) ,  denoted  as  A(f)  and  B(f) 
respectively,  are  orthogonal  linear  combinations  of  x(n),  it 
follows  that  A ( f )  and  B(f)  are  mutually  uncorrelated  Gaussian 
random  variables  each  having  the  distribution  N(0,  cA2)  .  The 
periodogram  of  x(n),  P,(f),  is  defined  as  the  sum  of  the 
squared  real  and  imaginary  parts  of  the  Fourier  Transform  of 
x (n)  t 

/’,(/)  =  a2(/)+ u3(/)  .  (2.ii) 


The  sum  of  the  squares  of  two  independent  zero-mean  normal 
variables  is  a  chi-squared  distribution  with  two  degrees  of 
freedom.  The  mean  and  variance  of  this  distribution  is  given 
by : 
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(2.12) 


4  oU 

P*  o,f 

8ctx4; 

r-  o.f 

(2.13) 


where  fN  denotes  the  sampling  frequency  [Ref.  4]. 

Proof  of  equation  2.13  for  frequencies  p^O,  N/2  is  given 
in  the  following  fashion: 

Consider: 

Px  =  A2  +  B 2 

where  A~N(o,crl)  (2.14) 

B~n\o,A 


We  know  that: 

(A2  +  B2f 
=  e[/44  +  2A2B2  +  B4] 

=  8al 


(2.15) 


(2.16) 


and  that  from  (2.12): 

E[p»)  =  2o«2 
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Therefore, 


Va-lr,)  =  E[/>,!]-(E[l'J,])J 
=  8ai-(2<4)J 
=  4oJ 


(2.17) 


F.  PERIODOGRAM  AVERAGING 

The  statistical  properties  of  the  periodogram  may  be 
improved  by  averaging  a  set  of  periodograms  together.  Assume 
that  K  independent  data  records  are  available,  all  for  the 
interval  o  <  n  <  (L  -  1)  and  all  are  realizations  of  the  same 
random  process.  The  data  is:  (x0(n),  0  <  n  <  L  -  1;  x,(n),  0  < 
n  <  L  -  1;  .  .  .  x,.,(n),  0<n<L-l}.  The  averaged 

periodogram  estimator  is  given  by: 

rAv{f)  =  TT  X1  'I'ERmif)  (2.18) 

where  I\er ,„(f)  t^le  Periodo9ram  of  the  mth  data  set: 


VrEn,„{f)  = 


L-i 

^r„,{n)exp{~j2nfnj 
n=0 


(2.19) 


The  mean  value  of  the  averaged  periodogram  will  be  the  same 
as  that  of  the  periodogram  based  upon  any  of  the  individual 
data  sets  since  periodograms  for  each  set  are  independent  and 
identically  distributed.  The  variance  of  the  periodogram  will 
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be  reduced  by  a  factor  of  K  as  a  result  of  the  averaging 
operation.  [Ref.  2] 

(2.20) 

In  actual  practice,  we  seldom  have  independent  data  sets.  It 
is  more  common  to  have  one  long  data  record  of  length  N.  A 
common  technique  is  to  segment  the  data  into  K  non-overlapping 
blocks  of  length  L,  where  N  =  KL.  Since  the  blocks  are 
contiguous,  they  cannot  be  uncorrelated  for  any  process  except 
white  noise.  Therefore,  the  actual  variance  reduction  is 
bounded  by  a  factor  less  than  or  equal  to  K.  If  the  data  are 
Gaussian  white  noise  samples,  the  autocorrelation  function  of 
the  data  will  decay  rapidly  and  the  blocks  will  be 
uncorrelated.  Thus,  the  periodograms  of  the  data  segments 
will  be  independent  and  (2.20)  will  be  accurate.  [Ref.  3]. 
As  an  illustration,  Figure  9  is  the  periodogram  of  64  samples 
of  Gaussian  white  noise  (zero  mean,  variance  1/2000).  Contrast 
this  with  Figure  10,  which  is  the  average  of  the  periodograms 
of  5  independent  64-point  data  records  obtained  by  segmenting 
a  320-point  record  of  white  noise  samples  with  the  same 
statistical  properties.  From  (2.11),  the  predicted  variance 
of  the  Figure  9  periodogram  is  4(l/2000)3  =  2.5  x  10"7  for  p  / 
0,  N/2.  From  (2.18)  we  would  expect  a  variance  reduction  by 
a  factor  of  1/N  =  1/5  or  6.9  dB  for  the  average  periodogram. 
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The  actual  measured  variance  reduction  between  the  single  and 
averaged  periodograms  is  6.7  dB.  A  variation  of  this 
averaging  scheme  was  proposed  by  Welch  [Ref.  5]  involving  the 
application  of  a  non-rect angular  window  function  to  each  data 
segment  and  overlapping  the  segments  (typically  in  a  4:1 
ratio) . 

In  interpreting  spectral  estimates,  it  is  important  to  be 
able  to  discriminate  between  spectral  detail  due  to 
statistical  fluctuation  and  actual  frequency  content.  A 
standard  way  of  evaluating  the  goodness  of  a  spectral 
estimator  is  via  confidence  intervals.  A  means  of  deriving 
a  confidence  interval  for  the  averaged  periodogram  is 
described  in  References  3  and  6. 
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G.  SPECTRAL  SMOOTHING  THE  DANIELL  PERIODOGRAM 

Daniell  suggested  that  a  means  of  smoothing  the 
fluctuations  of  the  periodogram  was  to  average  over  adjacent 
spectral  frequencies.  [Ref.  7]  He  proposed  a  modified 
periodogram  estimate,  P0(f),  in  which  each  frequency  spectral 
estimate  was  obtained  by  averaging  over  p  spectral  points  on 
both  sides  of  the  frequency  f  under  consideration.  The 
Daniell  Periodogram  is  given  by: 


1 

2p+l 


2  P{fn) 


w=i-p 


(2.21) 


A  generalization  of  this  concept  is  to  pass  the  sample 
spectrum  through  a  low-pass  filter  with  frequency  response 
H ( f ) .  The  Daniell  periodogram  may  then  be  expressed  as  the 
convolution  of  the  sample  spectrum  with  a  low-pass  filter  H(f) 
[Ref.  7]. 

I’oif)  ='’(/)*»(/)  ,2-22) 


The  larger  the  p  used,  the  greater  the  smoothing  effect 
will  be.  As  with  other  methods,  the  price  paid  for  smoothing 
is  a  loss  of  resolution.  Figure  11  shows  the  effect  of 
Daniell' s  operation  (p=2)  on  a  spectral  estimate  in  which  the 
frequency  of  the  test  signal,  10.0  Hz,  is  at  a  bin  center. 
Figure  12  shows  Daniell' s  method  performed  on  a  spectral 
estimate  where  the  frequency  of  the  test  signal,  10.7  Hz,  is 
not  at  a  bin  center. 
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In  summary,  the  FFT-based  spectral  estimation  technique 
(i.e.,  the  periodogram)  remains  popular  due  to  its 
computational  efficiency  and  good  performance.  Frequency 
resolution  (in  Hz)  is  proportional  to  the  reciprocal  of  the 
length  of  the  data  measured  in  seconds.  The  ability  to 
resolve  closely-spaced  signal  components  is  degraded  by  a 
combination  of  side  lobe  effects  and  main  lobe  broadening. 
Side  lobe  suppression  is  possible  through  the  use  of  non- 
uniformly  weighted  (non-rectangular)  window  functions  but  only 
at  the  cost  of  main  lobe  broadening.  Despite  these 
limitations  and  the  advent  of  modern  spectral  estimation 
techniques  such  as  parametric  modeling,  the  periodogram 
remains  the  most  popular  spectral  estimator  as  a  result  of  its 
relative  simplicity,  robustness,  and  ease  of  implementation 
in  real  time. 
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Figure  11.  Daniell  Spectral  Estimator,  p  =  2 
f  =  10.0  Hz,  64  points  Zero-padded  to  128 
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III.  KALMAN  FILTERING  IN  SPECTRAL  ESTIMATION 


A.  BACKGROUND 

A  continuing  problem  with  FFT-based  spectral  estimation 
schemes  is  the  trade-off  between  spectral  resolution  and  side 
lobe  suppression.  If  a  non-rectangular  window  function,  i.e., 
the  Hamming  or  Blackman  window,  is  applied  to  time  series  data 
for  the  purpose  of  minimizing  spectral  side  lobes,  the  side 
effect  is  a  loss  of  resolution  caused  by  the  broadened 
mainlobe.  In  general,  the  better  the  side  lobe  suppression, 
the  broader  the  main  lobe.  An  extreme  example  is  the  minimum 
4-sample  Blackman-Harris  window.  The  highest  side  lobe  of 
this  window  which  is  92  dB  down  from  the  main  lobe  peak.  The 
cost  of  this  level  of  side  lobe  attenuation  is  that  the  3-dB 
bandwidth  (main  lobe)  is  1.90  bins  versus  0.89  bins  for  a 
rectangular  window  [Ref.  1].  What  is  proposed  here  is  a  novel 
application  of  the  Kalman  filter  to  the  periodogram  for  the 
purpose  of  minimizing  spectral  sidelobe  effects  without  the 
usual  attendant  loss  of  resolution. 

The  Kalman  filter  program  demonstrated  here  was  written 
by  Dr.  Roberto  Cristi  at  the  Naval  Postgraduate  School, 
Monterey,  California  in  1988.  It  is  an  implementation  of  the 
filtering  algorithm  first  proposed  by  Kalman  and  Bucy  [Refs. 
8  and  9]  and  is  now  widely  used  in  control  system  theory. 
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Cristi's  program  was  originally  developed  to  detect  piecewise 
constant  segments  of  time  series  data  corrupted  by  noise. 

A  discrete  time  state-space  system  model  is  given  by: 


r(k  +  1)  =  *l>x(k)  +  Au«(fc)  +  A  pv{k) 


(3.1) 


y{k)  =  Cx(k)  +  w(k) 


where  x(k)  is  the  state  vector,  u(k)  is  the  input,  v(k)  is  an 
input  disturbance,  y(k)  is  the  observed  data  and  #(k)  is  the 
measurement  noise.  The  discrete  transition,  input,  input 
disturbance  and  observation  matrices  are  i,  A,,  Af,  and  C 
respectively.  The  input  disturbance  and  measurement  noise  are 
further  specified  by: 


£  [y(A;)y7  (k  +  w)j  =  Vri<5(m) 

£.[u;(/:)u;7  (/c  +  m)]  =  Wrf(5(,,i)  , 

where  and  W„  are  covariance  matrices. 

The  Kalman  gain  equations  are  given  by: 

L(k  f  Jj£)  =  H>r(k\kyi>T  +  AfVdAp 


(3.3) 


(3.4) 


(3.5) 


K(k  +  1)  =  F(k  -i  l|fc)C7' [l2P(k  +  \\k)Cr  +  Wrf]_I 


(3.6) 


r{k  -t  \\k  +  i)  =  [i  -  k(*  +  i)c]r(fc  1 1|*) 
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where  P(k+l|k)  denotes  the  covariance  matrix  at  time  k+i  given 
observations  to  time  k  and  K  is  the  Kalman  gain  matrix.  The 
Kalman  filter  equations  are  given  by: 

i(/c  +  l|/c)  =  <I>i(Jc)  +  A0u(/c)  (3.8) 

x[k  +  l|fc  +  l)  =  + 1|^)+  K(/c  +  l)Jy(fc  +  l)~Cr(fc  +  l|k)| 

where  x(k+l|k)  denotes  the  estimate  of  x  at  time  k+1  given 
observations  to  time  k.  Note  that  the  initial  condition 
P ( 0 | 0 )  must  be  specified  in  order  to  start  the  process: 

r(0|u)  =  t[(r(0)-i(0))(x(0)-i(0))r]  ^  (3.10) 

Equation  3.10  specifies  the  covariance  matrix  of  the 
initial  error.  The  covariance  matrix  is  a  measure  of  the 
confidence  on  the  initial  estimate  x(0) . 

Consider  the  simple,  one-dimensional  problem  of  detecting 
a  piecewise  constant  time  series  segment  corrupted  by  noise, 
which  was  the  original  purpose  of  Cristi's  program.  The 
signal  and  its  noisy  observation  are  given  by: 

x[k  +  1)  =  x(k) 
y(k)  =  x(k)  +  w(k) 

where  w(k)  is  the  corrupting  noise. 


(3.11) 

(3.12) 
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Now  define: 


'•(*)  =  #)-«(*)  (3.13) 
where  x(k)  is  the  estimate  of  x,  C  is  some  constant,  and  i  is 
the  innovation  sequence.  The  sequence  i(k)  represents  new 
information  not  contained  in  the  previous  observations  y(k- 
1)  ,  y (k-2) , . . .y(0) .  Elements  of  the  sequence  i  have  the 
property: 

*['(%(*  850  for  all  in  £  1  (3.14) 


Equation  3.14  states  that  each  element  of  i  is  orthogonal 
to  all  past  observations. 

Using  Baye's  theorem,  we  can  compute  the  probability  of 
the  observations  (y(k),  y (k-1) , . . .y ( 0) )  in  the  following 
fashion: 

\'t(y(k),y(k -  I)...y(0))  =  r. (y(*)|y(* -  1)...y(U))l'i(y(*  -  l)--y(O)) 

(3.15) 

Utilizing  the  recursive  property  of  this  expression,  we 
can  write: 


rr (y (k)> y ik  -  V- ■  ■  y (°))  =  FI  r*  (y  ( *)|y{  *  - 1 )•  •  ■ -y (o))  Vi (y (o));  k  z  1  ( 3 . 1  e ) 


Using  the  Orthogonality  Principle,  it  can  be  shown  [Ref. 

10]  : 

r,(y(%(* - J) -y(«)) - N(ci(k),a\k)cT + tv,)  ;  (3  17) 
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where  N  denotes  a  normal  distribution  and  W,  is  the  covariance 
matrix  of  the  observation  noise.  At  time  k,  it  is  then 
possible  to  compute  the  probability  Pr (y (k) | y (k-1) . . .y (0) ) . 
If  the  data  under  examination  consists  of  piecewise  constant 
segments,  then  at  each  new  observation  two  possibilities 
exist: 

1)  the  current  observation  is  a  continuation  of  the  last 
piecewise  constant  segment  of  data  observed  or 

2)  the  current  observation  is  the  first  element  of  a  new 
segment  of  data  with  a  constant  value. 

What  is  now  required  is  a  means  of  computing  the 

probability  that  a  transition  between  piecewise  constant 

sections  has  occurred.  Let  us  now  define  a  parameter  p  as 

a  means  of  quantifying  the  likelihood  of  a  transition  and  a 

binary  random  variable  1  as  follows.  If  a  transition  has  not 

occurred,  then  K  =  0  and  the  current  observation  is  filtered 

using  a  Kalman  filter  updated  with  the  current  gain.  If  a 
» 

transition  has  occurred,  then  ^  =  1  and  the  current 

observation  is  filtered  using  reinitialized  Kalman  filter. 
Now  define  the  probability  density  functions: 

I'’r(y(Jc)  =  0)  =  mexp(/J)  (3.18a) 

rr(r(k)  =  l)  =  mexp(-P)  ,  (3.18b) 
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_ _ 1 

where  m  ~  exp(/?)  +  exp(-/?)  and  k  denotes  the  time  index. 

Assume  that  each  'fc(k)  is  an  independent  event. 

N 

Pr(r(o),r(i)...r(M))=  XFr(r(k))  (3.19) 

K=0 


We  desire  to  maximize  the  expression  : 

Pr(y(^)|y(4^-i))  (3.20) 

where  y(k)  is  the  vector  of  observations  up  to  and  including 
the  current  time  k  and  jf_  (k-1)  is  the  vector  of  previous 
estimates  of  the  binary  random  variable  K  up  to  time  (k-1) . 
Equation  3.20  is  the  probability  of  a  transition  or  non¬ 
transition  (depending  on  =  0  or  2f  =  1)  ,  given  present  and 
previous  observations  y  and  estimates  of  .  We  desire  to 
maximize  Equation  3.20  with  respect  to  y(k)  and  If  (k-1)  where 

y(*)  =  [y(*),y(*-l)...y(0)] 

=  [y(4#-i)] 

and 

f(k  -  1  )  =  {f(k- 1),  7{k  -  2)...  ml  .  {3'22) 


39 


I©*> 


By  Bayes'  Theorem, 


Pr(r(*fe(fc)'£(*-1)) 

_  Pr[y(k)fy(k).y(k-1)) 
Pr(y(k),y(k-  1)) 

_  pr(y(fc),y(fc  - 1),  y(k),  y(k  - 1)) 
P*(y{k),Kk~l)) 


[pr(y(*)|y(fc  - 1),  y(k),  y(k  - 1) 

•Pr(j 

,(k-l),y{k),y(k-lj) 

Pr 

y(k),y(k-l 

» 

Assuming  that  #(k)  is  independent  of  y(k-l)  and  %  ( k— 1 ) , 
the  second  term  in  the  numerator  of  (3.22)  becomes: 

Pr(y(fc  - 1),  r(k),ftk  - 1))  =  Pr(r(k))l’r{y(k  -  i),r(k  - 1))  ( 3  ■ 24  ’ 

Equation  3.23  is  then  maximized  with  respect  to  y(k)  and 
(k-1)  by  the  expression: 

maxjpr 

‘  I. 

,  (3.25) 

=  maxjpr  (y(k) |#  - 1  ),Y{k),r{k  ~  l)p(y(*))j  f 
Define  the  likelihood  function: 

=  in[pr(r(fc)|#),y(fc-i))} 

=  ln{pr(y(Jt)|y(lt - 1), **).*<*  "  1))}  +  !"{*(**))}  .  (3'26> 
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Note  that  Pr (y (k) | y(k-l) ,  ?  (k)  ,  jf(k-l))  can  be  computed 
by  the  modified  version  of  (3.17): 

=r*(y(%(‘-i)<y(‘-2)-y(*-1)) 

~  N(Cx(k),Cl’(k)CT  +  Wd)  ( 3  • 2 7 ) 

where  2  is  the  time  interval  between  the  current  sample  k  and 
the  last  detected  transition.  Equation  3.27  is  evaluated  for 
the  two  cases  of  an  updated  or  reinitialized  Kalman  filter. 
The  probability  Pr('2'(k))  can  be  computed  via  (3.18).  Thus  it 
is  possible,  given  each  observation  and  those  proceeding  it, 
to  compute  the  probability  that  a  transition  between  constant 
valued  segments  has  or  has  not  occurred. 

By  selection  of  the  parameter  (3  (see  Equation  3.18),  it 
is  possible  to  adjust  the  likelihood  that  a  transition  will 
occur.  The  larger  the  (3  selected,  the  less  likely  the  filter 
is  to  reinitialize.  If  "too  small"  a  value  of  (3  is  selected, 
the  filter  will  reinitialize  too  often  and  little  smoothing 
of  the  data  will  be  done.  If  "too  large"  a  /3  is  used,  the 
filter  will  become  too  insensitive  to  fluctuations  in  the  data 
and  will  not  reinitialize  at  all.  In  this  case,  transition 
points  will  not  be  detected  and  the  original  data  will  be 
obliterated  (over-smoothed)  .  Thus  far,  (3  must  be  determined 
heuristically  depending  upon  the  type  of  data  under 
observation.  In  general,  noisier  data  (more  statistical 
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fluctuation)  will  require  more  smoothing  and  thus  larger 
values  for 

Figures  12  -  16  demonstrate  a  tost  of  the  Failuan  filter 
program  on  a  square  wave  of  amplitude  ±1  corrupted  by  Gaussian 
white  noise  of  variance  0.40.  Figure  13  shows  the  observed 
data  with  the  uncorrupted  signal.  Figures  14  -  16  show  the 
filtered  data  for  0  =  0.20,  4.00,  and  50.00.  Figure  14,  /?  = 
4.00,  shows  the  case  where  a  "good"  value  of  /3  has  been 
chosen.  Note  that  the  filter  correctly  detects  the  actual 
transitions  in  the  observed  data  and  reinitializes  only  at 
these  points.  As  a  result,  accurate  recovery  of  the  original 
waveform  is  achieved.  In  contrast,  Figure  15  shows  what 
occurs  when  too  small  a  /3  is  selected.  The  filter  becomes  too 
sensitive  to  noise  fluctuations,  mistakenly  interpreting  many 
of  them  as  transitions.  The  filter  reinitializes  too  often 
(see  lower  plot  of  transition  points)  and  less  than  optimum 
smoothing  is  performed.  Figure  16  is  the  case  where  too  large 
a  p  is  used,  rendering  the  filter  too  insensitive  to 
transitions  in  the  observations.  After  the  initialization, 
the  filter  never  detects  a  transition  and  thus  never 
reinitializes.  The  result  is  the  obliteration  (over 
smoothing)  of  the  true  waveform. 
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Figure  14.  Output  of  Kalman  Filter,  Square  Wave  (+1) 

Plus  Noise,  B  -  4.00 


B.  KALMAN  FILTERING  APPLIED  TO  THE  PERIODOGRAM 


Now  that  the  Kalman  filter  program  has  been  demonstrated 
on  a  simple  time  series,  the  question  arises:  Can  this 
algorithm  be  adapted  for  smoothing  spectral  data?  The 
objective  is  to  use  the  algorithm  to  smooth  the  periodogram 
spectral  estimate  with  minimal  broadening  of  the  main  lobe(s) 
of  the  dominant  spectral  responses.  Ideally,  an  appropriate 
value  for  the  parameter  0  is  selected  such  that  the  noise 
portions  of  the  periodogram  are  smoothed  and  transition  points 
are  detectable  on  either  side  of  the  spectral  main  lobe(s). 
The  end  result  is  a  smoothed  periodogram  with  the  narrow  main 
lobes  of  the  original,  unfiltered  periodogram  preserved.  The 
noise  "floor"  out  of  which  the  signal  peaks  rise  will  be 
better  defined  and,  hopefully,  the  frequency  resolution  of  the 
original,  unwindowed  periodogram  will  be  maintained. 

The  test  signal  used  is  a  single  sinusoid  (unit  amplitude) 
embedded  in  Gaussian  white  noise.  The  sinusoidal  frequency 
is  10.7  Hz,  which  is  not  at  a  bin  center.  The  signal  is 
sampled  at  64  Hz.  A  record  of  128  data  points  is  zero-padded 
to  256.  The  variance  of  the  additive  noise  is  varied  to 
create  input  (time  series)  signal-to-noise  ratios  (SNRs)  of 
-3,  -6,  -9,  and  -12  dB  where  SNR  is  as  defined  in  Chapter  II. 
Appendix  C  shows  10  different  noise  realizations  at  each  SNR 
for  a  given  value  of  0.  The  objectives  of  the  investigation 
were  three-fold: 
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1)  To  heuristically  determine  an  "optimum"  value  for  the 
parameter  p,  given  the  test  conditions,  at  the  different  input 
SNR ' s . 

2)  To  determine  the  input  SNR  of  the  time  series  (for  128 
delta  points  zero-padded  to  256)  at  which  the  Kalman  algorithm, 
given  the  "optimum"  p,  can  reliably  discriminate  noise 
perturbation  from  signal  peaks. 

3)  To  determine  if  the  Kalman  algorithm  preserves  the 
spectral  resolution  of  the  unfiltered  periodogram. 

After  many  trials,  it  was  determined  that  values  for  (3  in 
the  range  100,000  to  700,000  provided  the  best  compromise 
between  undersmoothing  and  oversmoothing  the  spectral  data. 
Within  this  optimum  range,  100,000  causes  the  least  smoothing 
and  700,000  the  most.  The  the  lowest  input  SNR  (time  series) 
at  which  reliable  signal  discrimination  was  achieved  was 
-6  dB.  At  -6  dB,  p  =  300,000  gave  generally  good  results. 
Signals  could  be  detected  at  SNRs  (time  series)  as  low  as 
-12  dB,  depending  on  the  noise  realization  (see  Appendix  C) . 

The  consequences  of  too  large  or  too  small  a  |3  in  the 
frequency  domain  are  analogous  to  the  time  series  example 
depicted  in  Figure  14  -  16.  Figure  17  illustrates  the  results 
of  the  Kalman  filter  at  an  input  SNR  (time  series)  of  -6  dB 
(128  data  points  zero-padded  to  256),  P=  300,000.  Note  that 
the  single  spectral  peak  due  to  the  sinusoid  has  been  left 
largely  unaltered  (unbroadened)  and  that  we  have  successfully 
smoothed  the  noise  portion  of  the  periodogram.  The  filtered 
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periodogram,  Figure  17,  more  closely  approximates  the  ideal 
model  of  a  spectral  peak  protruding  up  through  a  noise  floor 
of  constant  value.  In  all  cases,  the  Kalman  filter  was 
applied  to  periodograms  of  unwindowed  (rectangular  window) 
data.  This  resulted  in  the  most  narrow  of  possible  main  lobes 
and  provides  the  highest  resolution.  For  comparison,  a 
Hamming  window  was  applied  to  the  time  series  data  prior  tc 
transformation  (Figure  17)  .  Some  spectral  smoothing  is 
apparent  along  with  the  expected  main  lobe  broadening.  The 
noise  floor  is  far  less  apparent  than  in  the  Kalman  filtered 
periodogram.  Figures  18  through  20  demonstrate  the  effects 
of  varying  /3  for  a  given  noise  realization,  data/transform 
length  and  input  SNR.  In  Figure  18,  using  /?=  10.0,  we  obtain 
some  smoothing,  but  the  end  result  is  little  improvement  over 
what  is  obtained  with  the  Hamming  window  (Figure  18) .  Note 
that  even  at  this  low  value  of  f3 ,  we  have  smoothed  the  spectra 
and  preserved  the  narrowness  of  the  main  lobe.  Figure  19,  (3 
=  2.00  x  10',  illustrates  the  effect  of  a  /3  which  is  too  large 
for  the  given  input  SNR  and  noise  realization.  Note  the 
tapering  effect  on  the  higher  frequency  side  of  the  main  lobe. 
This  is  a  symptom  of  over-filtering  (over-smoothing)  caused 
by  too  large  a  value  of  0.  A  smaller,  closer-to-ideal  0  would 
have  caused  the  filter  to  reinitialize  after  the  peak  and  thus 
preserve  the  sharp  down-transition  of  the  original 
periodogram.  In  this  case,  the  filter  did  not  reinitialize 
and  smoothed  the  higher  frequency  side  of  the  main  lobe. 
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Whenever  this  tapering  effect  is  encountered,  better  results 
(sharper  main  lobe)  can  usually  be  obtained  by  reducing  p. 
Figure  20,  p  =  5.00  x  106,  demonstrates  obliteration  of  the 
original  spectra  caused  by  a  j3  which  is  grossly  too  large. 
Figures  B.l  and  B.2  in  Appendix  B  show  the  effect  of  varying 
p  over  a  wide  range  for  a  given  data  record  length,  transform 
length,  input  SNR,  and  noise  realization. 
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Figure  17.  Sinusoid  (f  «  10.7  Hz)  Plus  Noise,  SNR  -6  dn 

(a)  Periodogram  (rectangular  window) 

(b)  Periodogrom  (Hamming  window) 

(c)  Kalman  Filter  Output  p  =  300000 
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Figure  18.  Sinusoid  (f  =  10.7  Hz)  Plus  Noise,  SNR 

(a)  Periodogram  (rectangular  window) 

(b)  Periodogram  (Hamming  window) 

(c)  Kalman  Filter  Output  p  -  10.0 
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Figure  19.  Sinusoid  (f  -  10.7  Hz)  Plus  Noise,  SNR 

(a)  Periodogram  (rectangular  window) 

(b)  Periodogram  (Hamming  window) 

(c)  Kalman  Filter  Output  /3  =  2.00  E6 
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Figure  20.  Sinusoid  (f  -  10.7  Hz)  Plus  Noise 

(a)  Periodogram  (rectangular  window) 

(b)  Periodogram  (Hamming  window) 

(c)  Kalman  Filter  Output  (3  =  5.00  E6 


C.  EFFECTS  ON  SPECTRAL  RESOLUTION 


In  order  to  evaluate  the  effects  of  the  Kalman  filter  on 
spectral  resolution,  a  second  spectral  component  was  added  to 
the  test  data.  For  the  test  periodogram,  bin  width  is  fs/N  = 
64/128  =  0.5  Hz.  Note  that  we  used  N  =  128,  the  data  record 
size,  and  not  N  =  256,  the  transform  length.  As  stated  in 
Chapter  II,  zero-padding  does  not  improve  frequency 
resolution.  It  merely  allows  us  to  interpolate  more  frequency 
points.  Initially,  a  second  sinusoid  (also  unit  amplitude) 
at  13.9  Hz  was  introduced.  The  frequency  13.9  Hz,  like  10.7 
Hz,  is  not  a  bin  center  and  is  many  bin  widths  separate  from 
10.7  Hz.  With  /3=  30,000,  the  Kalman  filter  successfully 
discriminated  the  signal  peaks  from  the  background  noise  (see 
Figure  21) .  Next,  the  second  sinusoidal  frequency  was  brought 
in  to  11.2  Hz,  one  binwidth  separation  from  the  original 
signal  at  10.7  Hz.  This  is  close  to  the  0.89  binwidth 
resolution  limit  of  the  rectangular  window.  The  two  peaks  are 
clearly  visible  in  the  unfiltered  periodogram  (see  Figure  22)  . 
After  filtering  by  the  Kalman  filter,  the  spectral  estimate 
is  smoothed  and  the  resolution  of  the  original  periodogram  is 
preserved  as  evidenced  by  the  two  still-visible  spectral  peaks 
(see  Figure  22) . 
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D.  THE  NOISE-ONLY  AND  SIGNAL-ONLY  CASES 

The  effects  of  the  Kalman  filter  on  noise-only  and  signal- 
only  periodograms  was  tested.  Figures  23-25  show  the  Kalman 
filter  applied  to  three  different  realizations  of  Gaussian 
white  noise,  zero  mean,  0.5  variance.  As  before,  128  sample 
points  were  zero-padded  to  256.  Using  our  " ideal"  /3  of 
300,000,  no  sharp  spectral  peaks  were  discriminated.  This  was 
to  be  expected  since  no  dominant  spectral  component  was 
present.  Contrast  these  results  with  Figure  26,  which  is  the 
Kalman  filter  applied  to  signal-only  data.  In  Figure  26a,  the 
characteristic  sine  function,  translated  up  to  the  sinusoidal 
frequency  10.7  Hz,  is  visible.  Figure  26b  shows  the  well- 
known  smoothing  and  broadening  effects  of  the  Hamming  window. 
In  Figure  2Gc,  with  p  =  300000,  the  Kalman  filter  smoothed  the 
side  lobe  structure  of  the  sine  and  preserved  the  narrow  spike 
of  the  main  spectral  peak. 
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Figure  23.  Noise-Only,  Var  0.5,  Realization  1 

(a)  Periodogram  (rectangular  window) 

(b)  Periodogram  (Hamming  window) 

(c)  Kalman  Filter  Output  B  =  300000 
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Figure  24.  Noise-Only,  Var  0.5,  Realization  2 

(a)  Periodogram  (rectangular  window) 

(b)  Periodogram  (Hamming  window) 

(c)  Kalman  Filter  Output  (3  =  300000 
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gure  25.  Noise-Only,  Var  0.5,  Realization  3 

(a)  Periodogram  (rectangular  window) 

(b)  Periodogram  (Hamming  window) 

(c)  Kalman  Filter  Output  (3  =  300000 
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E.  THE  EFFECT  OF  DATA  RECORD  AND  TRANSFORM  LENGTH 

Finally,  the  number  of  data  points  was  increased  from  128 
to  256,  512  and  1024.  The  objective  was  to  evaluate  the 
performance  of  the  Kalman  filter  for  a  given  input  signal 
strength  (in  this  case  +12  dB)  at  different  length 
periodograms .  In  each  case,  the  data  record  was  zero-padded 
to  twice  its  original  length  (i.e.,  512  points  zero-padded  to 
1024)  .  Also  in  each  case,  the  input  SNR  was  decreased  in 
order  to  compensate  for  the  increased  processing  gain  caused 
by  the  data  record.  Processing  gain  is  approximated  by: 

processing  gain  [log,  (data  record  length)-l]x3  dB  (3.27) 

For  example,  in  our  baseline  case  of  128  points,  the 
expected  processing  gain  is  [log,  (128)-l]x3  dB  =  18  dB.  For 
an  input  SNR  of  -6  dB,  the  expected  output  SNR  is  then  18-6 
=  12  dB,  which  is  approximately  the  strength  of  the  peak  in 
Figure  17.  For  the  longer  data  trials,  the  additive  noise 
variance  (power)  was  increased  in  order  to  maintain  output  SNR 
at  approximately  12  dB.  Initial  results  indicate  a  dependence 
of  p  on  data/transform  length.  As  the  data/transform  length 
increases,  better  results  may  be  obtained  by  increasing  (3  (see 
Appendix  D) . 
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IV.  CONCLUSIONS 


The  Kalman  filter  can  enhance  the  spectral  peaks  of  a 
periodogram  of  an  unwindowed  time  series.  This  is  most 
apparent  in  the  single  spectral  peak  case.  In  the  case  of 
multiple  spectral  peaks,  the  resolution  of  the  unfiltered 
periodogram  is  largely  preserved  since  the  Kalman  filter  will 
smooth  the  spectral  estimate  without  major  broadening  the 
narrow  band  components.  Using  a  filter  parameter  in  the  range 
100,000  to  700,000  and  a  128-point  data  record  zero-padded  to 
256  points,  reliable  signal  detection  was  achieved  at  SNR's 
of  -6  dB  of  the  time  series.  Signal  detection  is  possible 
down  to  -12  dB  (of  the  time  series  SNR)  ,  depending  on  the 
noise  realization  used. 

Topics  for  further  study  are  the  application  of  the  Kalman 
filter  to  multidimensional  (time  varying)  spectra,  and 
quantification  of  selection  criteria  for  the  filter  parameter 
/3.  In  addition,  the  dependence  of  /?  o.i  input  SNR,  output  SNR, 
record  length  and/or  transform  length  should  be  examined. 
Another  possible  follow-on  project  is  the  development  of  an 
enhanced  Kalman  filtering  algorithm  that  adjusts  the  parameter 
P  based  on  the  assignment  of  signal  or  noise  only.  This  would 
mean  faster  filter  response  during  signal  portions  and  slower 
response  during  noise-only  segments. 
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APPENDIX  A 
COMPUTER  CODE 

The  Kalman  filtering  program  was  originally  written  in 
FORTRAN  77.  The  FORTRAN  code  is  given  in  Appendix  A.l.  For 
this  thesis,  the  filter  program  was  converted  to  PC-MATLAB 
(Version  3.13)  and  simulations  run  on  an  80386-based  IBM 
compatible  PC.  The  MATLAB  code  for  the  filtering  program  is 
given  in  Appendix  A. 2 


SECTION  A. 1 
FORTRAN  Computer  Code 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

5  0 
2  2? 


c: 


66  6 


c 


**  Nonstationary  filtering  using 
**  subopt imal  kalman  filtering  (local 
**  average),  and  gibbs  field  with 
**  anihiling. 

**  The  input  file  must  given  on  INPUT . DAT 
**  The  filtered  output  file  is  stored  on  OUTPUT.DAT 
**  The  detected  breackpoints  are  given  in  MODEL. DAT 
**  All  these  * . DAT  files  are  ASCII. 


**  The  program  now  works  for  12B  data  points  (see  the  variable 
**  "npoints"  below.  This  can  be  changed  to  any  number  of  points. 


*  * 
*  * 
«  * 
*  * 
*  * 
4  * 
4  * 
4  4 
4  4 


The  program  requires  to  enter  2  parameters: 

Msigma**:  the  value  of  the  noise  standard  deviation  (nonzero); 

"beta  •*:  a  positive  parameter.  It  is  a  measure  of  the  probability 
the  signal  having  a  jump.  As  is  now  this  parameter  is 
set  by  pure  trial  and  error.  If  you  get  too  many 
jumps  detected  it  means  that  beta  is  too  low.  If  you  get 
too  few  jumps  it  means  that  beta  is  too  large.  However 
usually  the  best  value  of  beta  depends  on  the  signal  to 
noise  ratio  of  the  data. 


real  yin(256),  y(256),  x(2,256) 

real  kl,k2 

integer  pointer(2,  256),  t,  out 
integer  mout(256) 

open(l,  f i 1 e= ' output . dat ' ,  st a tus= ' ol d ' ) 

open  (2,  file=  # input .dat # ,  status=/old ' ) 

open(3,  file=  'model. dat',  s t a t us= 'old ' ) 

**  get  data  from  file 

rewind  1 

rewind  2 

rewind  3 

4  4 

npo 1 nt s=l 28 

4  4 

do  50  t-l,npoints 

read  (2,222)  y(t) 

yin(t)»y(t) 

wr i te ( * ,  ill)  y(t) 

con  t  i  ntie 

format ( f  d . 4 ) 

**  enter  data  and  initialize 
write(*,555) 

format ('  ENTER:  sigma, bet  a  (>o)') 
read(*,666)  sigma,  beta 
format ( 2  f 1 0 . 4 ) 
sv2  =  sigma  *  *  2 

e 1 "0 . 0 
e2=0 . 0 
d 1 -0 . 0 
d  2  =  0 .  o 
x  ( i  / 1 ) =y  ( i ) 
x  ( 2 , 1 ) -y ( 1 ) 
t  a  u 1 


*  *  main  1 oop 

do  100  t- 1 ,  (npoints-1) 
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kl«l.O/ (tau 41 . 0) 

xll«x(l# t)  +  kl* (y ( t) -x( 1 , t ) ) 

dl l=dl-beta 

ell«el+ (1. 0/ (2 ,0*sv2) ) * (y (t+1) -xll)  **2 
cl  l«*el  l4dl  1 
c 

k2®1.0 

xl  2=x ( 1 , t )  4  k2  * (y ( t) -x ( 1 #  t ) ) 
d 1 2=d 1 4  beta 

el2=el4(1.0/(2.0*sv2))*(y<t4l)-xl2)**2 

Cl2=el24dl2 

c 

k  1  -  0 . 5 

x2 1 =x ( 2 , t )  4  kl * (y ( t i -x ( 2  f  t ) ) 
d2  1  =d2 -be t a 

e21=e2  4 ( 1 . 0/ ( 2 . 0*sv2 ) ) * ( y ( 1 4 l ) -x2 1 ) *  * 2 
C21=e2l4d21 
c 

k  2  =  1 . 0 

x22  =  x ( 2  ,  t )  4  k2*(y(t)-x(2,t)) 
d2  2-d2  4  bet  a 

e22-e24 (1.0/(2.0*sv2) ) * (y(t*l)-x22)  *  *  2 
c2  2  =  e2  2  4  d2  2 
c 

c  wr i te ( * , 444 )  t , d 1 1 , el 1 , d 1 ? , el  2 , d2 1 , e2 1 

4  44  format (i3,3(fl0.2, f 10.2,2x) ) 

c 

c:  ***  update  states  in  dynamic  prog, 

i f (cl  1 . 1 e . c2 1 )  then 
x ( 1 , t » 1 ) =  xl  1 
el=el 1 
d  1  -  d  1  1 
*  cl =o 1 4d 1 

pointer  ( 1 ,  t  ♦  1 )  =1 
t au  =  tau  4 l 

else 

X  ( 1 , t  4 1 ) -x2  1 

e  1  -  e  2  1 
dl =d2 1 
cl  =  e 1 4  d 1 
tau  ~ 2 

poi  liter  ( 1 , 1 4  l )  =2 

end  i  f 

c 

if  (c22.lt.cl2)  then 
x(2,tU)*x22 
e2=e22 
(12-122 
c  2  -  °  2  *  d  2 

pointer(?;t^i>'? 

else 

x  (  2 , t  4 1 ) -x  1  2 

e2  e) 2  : 

d  2  -  d  1  2 
c  2  e  2  4  d  2 
point'u  (?  ,t  '  1  )  :  1 


6  / 


end  i  f 


continue 


100 

C 

c 


c 


e 


c 

c 

333 

150 

C 


3  34 


0  00 


777 
1  1  1 
c 


backward  substitution  and  smoothing 
tau=»l .  0 

i f (Cl . 1 e . c2 )  then 
out  =  l 


else 
end  i  f 


out=2 


y (npoints) =  x (out  # npoints) 
n  2  =  0 

do  150  t=npoints, 2, -1 

out=pointer(out,t) 
xout=x (out , t-1 ) 
if  (out. eg. 2)  then 
tau  =  l . 0 


else 

tau=tau *1.0 

end  i  f 


y (t-l)*y (t) ♦ (1.0/tau) * (xout-y (t) ) 

mout ( t - 1 ) =out  *  1 0  0 

write  ( 1 , 1 11 )  xout 

n2=n2 ♦ out-1 

wr i te ( * , 333)  t,  out 

format ( 2 ( 2x , i 5) ) 

continue 


s 1 gma=0 . 0 

do  800  t=l, npoints 

write(l,  HI)  y(t) 

write(3,334)  mout(t) 

formats  i  5 ) 

ye~(y (t) -yin(t))**2 

sigma=sigma  t  ( 1 . 0/t ) * (yo-s i gma ) 

conti nue 

sigma=sqrt (sigma) 

write(*,777)  sigma,  n2 ,  npoints 

format ( 9  sigma-',  f8.4,  '  n2=',I5,'  npoi nts- ',15) 

format (  f  8 . 4 ) 

rewind  1 
rewind  2 
stop 
end 


SECTION  A. 2 

PC-MATLAB  Computer  Code 


EC  THESIS  GO,  W.  W. 

THE10.M  KALMAN  FILTER  APPLIED  TO  PERIODOGRAM  OF  TWO  SINUSOIDS 
IN  GAUSSIAN  WHITE  NOISE.  128  DATA  POINTS  ARE  ZERO- 
PADDED  TO  256  AND  THEN  THE  PERIODOGRAM  IS  COMPUTED. 

ONLY  HALF  OF  THE  RESULTING  FREQUENCY  POINTS  (UP  TO 
ONE-HALF  OF  THE  SAMPLING  FREQUENCY)  ARE  PLOTTED  AND  USED 
AS  INPUT  TO  THE  KALMAN  FILTER.  THE  FOLLOWING  CASES 
ARE  PLOTTED: 

1)  PERIODOGRAM,  RECTANGULAR  WINDOW  ON  TIME  DATA 

2)  PERIODOCRAM,  HAMMING  WINDOW  OH  TIME  DATA 

3)  OUTPUT  OF  KALMAN  FILTER  APPLIED  TO  PERIODOGRAM 
OF  RECTANGULARLY  WINDOWED  DATA  (CASE  1). 


The  program  requires  2  parameters  to  be  specified: 

"sigma":  the  value  of  the  noise  standard  deviation  (nonzero); 

"beta  ":  A  positive  parameter.  It  is  a  measure  of  the  probability 
of  the  signal  having  a  jump.  Now  this  parameter  is 
set  by  pure  trial  and  error.  If  you  get  too  many 
jumps  detected  it  means  that  beta  is  too  low.  If  you  get 
too  few  jumps  it  means  that  beta  is  too  large. 

NOTE  1:  THIS  PROGRAM  UTILIZES  MATIJVB  FUNCTIONS  PER .  M  AND  PERLN.M 
(CODE  FOLI.OWS  MAIN  PROGRAM)  TO  COMPUTE  THE  PERIODOGRAM 
IN  dB  AND  LINEAR  UNITS  RESPECTIVELY. 

NOTE  2:  THIS  FROGPAM  UTILIZES  MATIAB  FUNCTION  FVEC.M  (CODE 
FOLIOWS  MAIN  PROGRAM)  TO  CREATE  A  FREQUENCY  VECTOR 
FOR  PLOTTING . 


fl  =  10.7Hz,  MOT  A  BIN  CENTER 
f 2-  11.2,  NOT  A  BIN  CENTER 


fs  =  6 4  Hz,  SAMPLING  FREQUENCY 
128  DATA  POINTS  ZERO-PADDED  TO  256 
WHITE  NOISE  VARIANCE  =  4000/2000 
INPUT  SNR  -6.02dH 


clear 

clg 

f  1  =  10.7  ,' 

f  2  -•  11.2 
fs  6  4  ; 


I 

%  f  is  frequency 
1  fs  is  sampling  frequency 


nvar=400 0/2000; 


1  noise  variance 


for  n=  0  :  127  ;  l  compute  signal  vector 

x(n+l)  =  cos ( n * 2 *p i * ( f 1 / f s) )  +  cos ( n* 2 *pi * ( f 2/ f s ) )  ; 

end 


i 

rand ( ' norma  1  '  )  ; 
t  and ( ' seed ' , J  ) 

nz=sqrt ( nvar )  . *rand ( 1 : 1 20 )  ;  1  noise  vector 

xn~x  +nz;  1  corrupt  signal  with  noise 
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w=hammlng ( 128 )  ; 
xnw-w'  .  *xn; 


I  Hamming  window 
%  apply  Hamming  window 


xp=  [  xn  zeros(l:128)  ); 
xpw=  {  xnw  zeros(l:128)  )  ; 

psd  -perln(xp);  %  periodogram(linear  units) 

test-per (xp) ;  t  per iodoyrar ;dB) 

testw=per (xpw) ; 

f req=fvec(64 , xp) ;  \  frequency  vector  for  plotting 

subplot (211) ,plot(freq(l:128) ,test(l:128) ) 

t i tie ( 'THE10 : 2  SIN  IN  NOISE, SNR  -6 . 02dB, f 1=10 . 7 , f 2=1 1 . 2 ' ) 
xlabel (' frequency ' ) 
ylabel ( 'magnitude' ) 

subplot (211) ,plot(freq(l:128) , testw( 1 : 128 ) ) 

t  i  1 1  e  (  '  THE  1 0  :  2  SIN  IN  NOISE,  HAM  WIN, SNR  -  6 . 0  2dB ,  f  1  =  1  0 . 7  ,  f  2  -  1 1  .  2  '  ) 

xlabel ( ' frequency ' ) 

ylabel ( ' magnitude ' ) 

meta  preplt2 

pause 


y=  psd (1:128)  ; 


KAI.MAN  FILTER 

y  is  DATA  RECORD.  FILTER  IS  AFTLIED  TO  TERIODOGRAM  IN 
LINEAR  UNITS. 

x=zeros (2,128) ; 
pointer-zeros(2,  128)  : 
yin  =  ys 

beta  -500000.0; 
sigma  *=  sgrt(nvar); 
sv2  =  s igma *  2 ; 

npoints=length (y) ; 

el-0 . 0; 
e  2  =  0 . 0  ; 

d  1  =  0 . 0  ; 
d  2  -  0 . 0  ; 

x ( l , i ) -y ( l )  ; 

X ( 2 , 1 ) =y  ( 1 )  ; 

tau  =  l . 0 ; 

MAIN  LOOP 


%  filter  parameter. 

I  noise  standard  deviation 
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for  t-1 : (npolnts-1)  ; 

kl=l . 0/ ( taui 1.0); 
xll=x(l,t)+kl*(y(t)-x(l,t)); 
dl 1 =dl-beta ; 

ell=el+ (1.0/  (2 . 0*sv2) ) * ( (y ( t+ 1 ) -xll ) " 2 )  ; 
cll=ell+dll ; 


k2=l . 0 ; 

xl2=x(l,t)<k2*(y(t)-x(l,t)); 
dl2=dl+beta ; 

el2=el+ ( 1 . 0/ (2 . 0»sv2 ) ) * ( (y ( t+1 ) -Xl2 ) ' 2 ) ; 
Cl2-»el2  +  dl2  ; 


kl=0. 5; 

x21=x(2,t)+kl*(y{t)-x(2, t) ) ; 
d21=d2-beta: 

e2 1  -  e2i (1.0/(2.0*sv2))*((y(t+l)-x21)*2) 
c21=e21+d21 ; 


k2-l . 0 ; 

x22=x(2,t)<k2*(y(t)-x(2,t)); 
d22=d2r  beta : 

e22=e21(1.0/(2.0*sv2))*((y(t+l)-x22)*2); 
C22=e22*d22 ; 

%  UPDATE  STATES  IN  DYtlAMIC  PROGRAM. 

if  cl  1 <c2 1 

x(l,ttl)-xll; 
el=ell  ; 
dl-dll  ; 
c  1  -  c  1  t  d  1  ; 

poi liter  ( 1 ,  t  *  3  )  =1  ; 

tau  =  taui 1 ; 

else 

x ( 1 . t i 1 ) =x2 1  ; 

e 1 =e2 1 ; 
d  1  -  d  2  1  ; 
cl ~el i dl ; 
tan- 2  ; 

po  i  liter  ( l  ,  t-t  1 )  =2  ; 

end 


if  c22<ci: 

x { ? , t  *  1 ) =x2  2  ; 
e  2  --  e  2  2  ; 
d  2  -  d  2  2  ; 
c2^e2 id?  ; 
pointer (2 ,ti 1 )  -  2  ; 

else  i 

y (2 , 1 1 1 ) -xl2  ; 
e2-el 2  : 
d  2  -  d  1  2  ; 
c  2  -  e  2  i  d  2  ; 
pointer ( 2 , t M ) =1 : 
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end 


end 

END  MAIN  LOOP 

BACKWARDS  SMOOTHING  AND  SUBSTITUTION 


tau=l . 0 ; 
if  cl<c2 

out=l ; 

else 


out=2 ; 


end 

y (npoints) =x (out , npoints) ; 

for  t»128:-l:2 

out=pointer(out, t) ; 
xout=x (out , t-1 ) : 


if  out==2 

tau=l . 0 ; 


else 


y ( t-1 ) =xout ; 
tau=tau  *  1  ; 


end 


y(t-l)=y(t)i (1.0/tau) * (xout-y ( t ) ) : 

y ( t-1 )  =xout; 

end 


trans ( t-1 ) =out  r 

end 

ynorm=(l/max(y) ) .*y: 
ydb=10* loglo (ynorm) ; 

ysh=  (  ydb (2 : length ( ydb) )  ydb(i)  ); 

subplot (212) .plot ( freq( 1 : 128) , ysh) 
t i t le ( 'THE 10 : KAL, SNR  -6.02, BETA  500000.0’) 

xlabcl ('frequency') 
ylabel ( 'magnitude  dB') 

meta  prepltl 
%  pause 

i  plot (trans, '  +  '), tit le (' transit  ion  pts') 
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%  EC  THESIS  GO.  W.W. 

%  PER.H  COMPUTE  THE  PERIODOGRAM  OF  DATA  VECTOR  X 


function  y*per(x) 
l=length (x) ; 
tr=f ft (x) ; 

for  1=0: (1-1)  ; 

ps(i  +  l)  =  (abs(tr(i-U) ) )  *2; 

end 

psnorm= ( 1/max (ps) ) . *ps ; 
y=10*log!0 (psnorm) ; 


EC  THESIS 


GO.  W.W. 


t  PERLN.M  COMPUTE  THE  PERIOUOGRAM  OF  DATA  VECTOR  X 

\  LINEAR  UNITS 

function  y=perln(x) 
l=length(x) ; 
tr=f f t (x) ; 

for  i=0: (1-1) ; 

ps(i+l)-=(abs(tr(l4))))"2; 

end 

y=ps; 


7-1 


\  EC  THESIS 

\  FVEC.M 

\ 

\ 

function  f 


GO , W . w . 

CREATE  THE  FREQUENCY  VECTOR  USED  IN  PLOTTING 
A  FERIODOGRAM.  fs  IS  THE  SAMPLING  FREQUENCY 
AND  X  IS  THE  DATA  VECTOR. 

=f vec ( fs , x) 


n= length (x)  ; 
£=fs* (0:n-l)/n; 


APPENDIX  B 

EFFECTS  OF  THE  KALMAN  FILTER  PARAMETER 

The  effects  of  changing  the  parameter  P  on  the  performance 
of  the  Kalman  filter  were  investigated.  The  test  data  was  a 
single  sinusoid,  with  frequency  of  10.7  Hz,  embedded  in 
Gaussian  white  noise.  The  frequency  10.7  Hz  was  specifically 
chosen  so  as  not  to  be  at  a  bin  center  of  the  FFT.  A  record 
of  128  data  points  was  zero-padded  to  256.  The  input  SNR  of 
the  time  series  was  -6  dB.  The  same  noise  realization  was 
used  for  all  runs.  Figure  B.l  shows  the  unfiltered 
(rectangular  windowed)  and  Hamming  windowed  periodograms . 
Figures  B.2a  through  B.2j  show  the  filtered  periodograms  for 
ten  different  values  of  (3.  As  discussed  in  Chapter  3,  p  in 
the  range  100,000  to  700,000  produced  the  best  results. 
Values  for  (3  below  this  range  tend  not  to  smooth  the  spectral 
estimate  enough  to  signif ic<antly  enhance  the  main  spectral 
peaks.  Values  of  p  above  this  range  tend  to  oversmooth  and 
obliterate  the  spectral  estimate  (depending  on  the  noise 
realization) . 


magnitude 


frequency  Hz. 
(a)  Rectangular  window 


(b)  I’eiiodogram  (Hamming  window) 


Figure  B.l.  Unfiltered  Feriodograms 


magnitude  dB  •  magnitude  dB 


frequency  Hz 


(a)  p  =  10.0 


10  20  30  40 

frequency  f  Iz 


(b)  P  =  100.0 

Figure  B.2.  Output  of  the  Kalman  Filter 


(c)  p=  1,000.0 


KAL  OUT, SNR  -fitlU, HIITA  10000.0 

0  f- - - r  t - 1 


!rc(|uciu\  Hz 


(d)  p  =  10,000.0 


Figure  B . 2 . 


output  of  the  Kalman  Filter  cont . 


mannitude  dB  magnitude  dB 


0 


-10 


KAL  OUT, SNR  -6dB,BF:TA  3(KK)00.0 


'  -20 
-30 


0  10  20  30 

liei|iiency  llz 


40 


(f)  0  -  300,000.0 


Figure  B.2.  Output  of  the  Kalman  Filter  cont 
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magnitude  dB  magnitude  dB 


frequency  I  fz 


(g)  /3=  700,000.0 


frequency  Hz 
(h)  0  =  1.00  E6 


Figure  B.2.  Output  of  the  Kalman  Filter  cont 


magnitude  dB  magnitude  dB 


frequency  Hz 


(i)  p  =  2.00  E6 


frequency  Hz 


(j)  p  =  5.00  E6 


Figure  B.2.  Output  of  the  Kalman  Filter  cont 


APPENDIX  C 

PERFORMANCE  OF  THE  KALMAN  FILTER  AT  DIFFERENT 
INPUT  SNR'S  ON  MULTIPLE  NOISE  REALIZATIONS 

The  performance  of  the  Kalman  filter  at  different  input 
SNRs  was  evaluated.  The  test  case  was  a  single  sinusoid, 
frequency  10.7  Hz  (not  a  FFT  bin  center).  A  record  of  128 
data  points  was  zero-padded  to  256.  The  Kalman  filter  was 
run  on  data  with  time  series  SNRs  of  -3,  -6,  -9,  and  -12  dB. 
Ten  different  noise  realizations  were  used  at  each  SNR.  Plots 
are  shown  in  Figures  C.l  through  C.40.  For  comparison,  the 
unfiltered  and  Hamming  windowed  pe_iodograms  are  also  shown 
for  each  simulation.  At  -6  dB  (time  series  SNR) ,  reliable 
detection  was  achieved  for  all  noise  realizations  tested. 
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REC  WIN.SNR  -3dB.f=  10.7  »  HAM  WIN.SNR  -3dB,f- 10.7 


frequency  Hz 


REC  WIN, SNR  -3dB,f=  10.7 
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ppniiuSinu 
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Figure  C.3 


-3  dB  Input  SNR,  Noise  Realization  3 


aprujimmu 


.■jpniiuSiHu  £jp  ppmuiriiHu 


Figure  C.4.  -3  dB  Input  SNR,  Noise  Realization  4 
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HAM  WIN, SNR  -3dB,f=10.7 


opniuirsnui 


f]p  spniiuiiniu 


Figure  C.7.  -3  dB  Input  SNR,  Noise  Realization  7 


REC  WIN.SNR  -3dB,f=10.7  n  HAM  WIN, SNR  -3dB, 


Figure  C.9.  -3  dB  Input  SNR,  Noise  Realization  9 
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-3dB. 
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Figure  C.12.  -6  dB  Input  SNR,  Noise  Realization  2 
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ncv 


frequency  frequenc 


freti  ue  ncv 


KAL  OUT.SNR  -6.BETA  300000.0 


cin  ppniiuSiHU 


frequency 


N.SNR  -6dB.f=  10.7  .  HAM  WIN.SNR  -6dB,f=10.7 


,SNR  -6dB, 


RE C  WIN. SNR  -6dB.f  =  10.7  „  HAM  WIN.SNR  -6dB,f=10.7 


-6  dB  Input  SNR,  Noise  Realization  8 
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Figure  C. 18 


frequency 


REC  WIN. SNR  -OdB.f-  10.7  „  HAM  WIN.SNR  -6dB.f=  10.7 
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-6  dB  Input  SNR,  Noise  Realization 


Figure  C . 19 


frequency 


-9  dB  Input  SNR, 
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Figure  C.20 


Noise  Realization  10 


Irequencv 


HAM  WIN, SNR  -9dB,f=10.7 


spnjiuSiuu  gp  opnjiuSiuu 


-9  dB  Input  SNR,  Noise  Realization  4 
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Figure  C.24 


frequency  f  fz 


apnjiuSeui 


apruiuSciu 


gp  apnjruthjui 


Figure  C.25.  -9  dB  Input  SNR,  Noise  Realization  5 
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Figure  C.27. 


-9  dB  Input  Sf 
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frequency  frequency 


gp  spnjiuScm 


Noise  Realization  7 
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Figure  C.29.  -9  dB  Input  SNR,  Noise  Realization  9 


112 


REC 


opnjjuSiuu 


gp  apnjuiSiHU 


Figure  C.30.  -9  dB  Input  SNR,  Noise  Realization  10 
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re  c]  ue 


ppruiunniu 
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Figure  C.32. 


-12  dB  Input  SNR,  Noise  Realization  2 


Figure  C.35.  -12  dB  Input  SNR,  Noise  Realization  5 
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RRC  WIN.SNR  -12dB.f=  10.7  n  HAM  WIN, SNR  -12dB,f=10.7 


spnjiuSuiu 
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Figure  C.36. 


-12  dB  Input  SNR,  Noise  Realization  6 


REC  WIN,SNR  -12dB.f=  10.7  „  HAM  WIN.SNR -12dB,f=  10.7 
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Figure  C.38.  -12  dB  Input  SNR,  Noise  Realization  8 
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spruiufjcui  gp  apnjiuSmu 


Figure  C.39.  -12  dB  Input.  SNR,  Noise  Realization  9 
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REC  WIN.SNR  -12dB,f=  10.7  „  HAM  WIN, SNR  -12dB,f=10.7 


Figure  C.40.  -12  dB  Input  SNR,  Noise  Realization  10 
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APPENDIX  D 

EFFECTS  OF  DATA/TRANSFORM  LENGTH 

The  effects  of  varying  the  data/ transform  length  on  the 
performance  of  the  Kalman  filter  were  investigated.  The 
baseline  test  case  was  a  pair  of  sinusoids,  frequencies  10.7 
and  13.9  Hz  (not  at  FFT  bin  centers)  embedded  in  Gaussian 
white  noise.  Data  records  of  128,  512,  anc  1024  points  were 
zero-padded  to  twice  their  original  length.  As  discussed  in 
Chapter  III,  the  time  series  input  SNR  was  decreased  with 
increasing  data/transform  length  in  order  to  compensate  for 
the  higher  processing  gains  of  the  longer  data  records  (so  as 
to  maintain  the  SNR  at  the  input  to  the  Kalman  filter  at 
approximately  12  dB)  .  A  /?  of  300,000  was  used  since  this 
value  gave  good  results  with  the  baseline  128  data  point  test 
case.  Results  are  given  in  Figures  D.l  through  D.3.  For 
comparison,  the  unfiltered  periodograms  are  also  shown  for 
each  data/transform  length.  Results  indicate  that  as  data 
transform  length  is  increased,  (3  may  have  to  be  increased  in 
order  to  obtain  optimum  smoothing  of  the  spectral  estimate. 
The  dependence  of  /3  upon  transform/data  length  is  a  potential 
topic  for  follow-on  study. 
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RFC  WIN. 2  SIN  IN  NOISF.,512  DATA  PTS,SNR  -12dB,tT  =  10.7,f2=13.9 
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Figure  D.2. 

512  Data  Points 
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f2=  13.9 
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